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OPTIMAL  HANKEL-NORM  APPROXIMATION 
AND  RATIONAL  FUNCTION  MODELS  USING 
IMPULSE  RESPONSE  DATA 


INTRODUCTION 


It  is  well  known  that  the  underlying  mathematical  formulations  in  systems  theory  and 
signal  processing  are  quite  similar  [1,2].  In  particular,  the  transfer  function  H  ( ~ }  with  j 

the  sequence  of  unit  impulse  responses  {/in},  plays  the  central  role  of  these  two  important  ;  '  \ 

areas  of  research.  In  the  ideal  situation  when  the  sequence  { hn }  is  noise-free,  then  it  is  very  \  ?  -v  ( 

easy  to  obtain  the  rational  function  model  H(z).  One  method  is  the  use  of  Pade  approxi-  _ - 

mants.  This  method,  being  linear,  is  particularly  simple.  In  systems  theory,  however,  the 
dimension  of  the  system,  which  agrees  with  the  degree  of  the  denominator  of  the  rationed 
function  H(z),  although  finite,  may  be  unreasonably  high.  Hence,  many  model  reduction 
methods  have  been  introduced  in  the  literature.  All  these  methods  amount  to  approxima¬ 
tion  of  the  rational  function  H(z)  by  another  rational  function  with  a  much  lower  •  egree 
(or  lower  system  dimension).  However,  only  one  method  stands  out  to  be  the  most  desir¬ 
able.  This  method,  known  as  optimal  Hankel-norm  approximation,  has  three  outstanding 
features:  first,  the  lower-degree  rational  function  approximant  is  uniquely  determined  by 
an  optimal  criterion;  second,  this  rational  function  approximant  is  guaranteed  to  be  sta¬ 
ble;  and  third,  the  exact  measurement  of  the  error  of  approximation  can  be  determined. 

The  mathematical  description  and  derivation  of  the  optimal  Hankel-norm  approximation 
is  especially  intriguing.  Based  on  the  infinite-dimensional  theory  due  to  Adamjan,  Arov, 
and  Krein  [3],  better  known  as  AAK,  it  relates  best  rational  approximation  in  the  supre- 
mum  norm  and  optimal  approximation  of  the  corresponding  Hankel  operator  by  Hankel 
operators  with  specified  finite  rank.  The  intimate  relation  between  rational  functions  and 
finite-rank  Hankel  operators  is  governed  by  the  beautiful  classical  result  due  to  Kronecker 
(cf.  Ref.  4);  and  the  relation  between  best  rational  approximation  in  the  supremum  norm 
and  optimal  Hankel  approximation  in  the  operator  norm  is  initiated  by  a  fundamental  re¬ 
sult  due  to  Nehari  [5].  The  AAK  theory  may  be  considered  as  a  generalization  of  Nehari’s 
theorem  in  that  the  singular  values  are  used  to  give  the  exact  measurement  of  the  error  of 
approximation.  Note  that  the  first  singular  value  agrees  with  the  spectral  radius,  which  is 
the  same  as  the  distance  in  the  supremum  norm  of  the  transfer  function  from  the  Hardy 
space  H°°  of  bounded  analytic  functions  in  the  open  unit  disc.  This  is  the  main  theorem 
due  to  Nehari. 

Unfortunately,  the  AAK  theory  is  very  deep,  and  there  seems  to  be  no  easy  way 
to  find  the  AAK  optimal  approximants.  Nevertheless,  when  the  transfer  function  is  al¬ 
ready  a  rational  function,  S.Y.  Kung  [6]  gave  an  algorithm  to  compute  the  AAK  optimal 
approximation.  This  method  is  called  control  (or  systems  reduction  via  Hankel  ap¬ 
proximation)  in  systems  theory.  However,  even  in  the  absence  of  noise,  impulse  response 
data  from  underwater  acoustic  transducers  do  not  yield  a  rational  model;  and  in  general, 
data  obtained  from  a  rational  model  are  contaminated  with  noise.  Hence,  for  all  practical 
purposes,  the  transfer  function  to  be  identified  is  not  in  the  form  of  a  rational  function,  and 
a  rational  approximation  criterion  is  necessary.  Recently,  Chui,  Li,  and  Ward  [7,8]  proved 
that  if  the  truncated  Hankel  operators  are  used  in  place  of  the  original  Hankel  operator 
(which  represents  the  transfer  function  to  be  identified),  then  the  AAK  approximants  of 
the  truncated  operators  converge  to  the  AAK  approximant  of  the  desired  transfer  func¬ 
tion.  We  will  derive  a  computational  procedure  based  on  this  convergent  result.  In  other 
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words,  if  a  sequence  {/in}  of  unit  impulse  responses  is  given,  we  will  derive  a  computational 
scheme  to  find  the  stable  rational  models  that  are  optimal  in  the  sense  of  best  Hankel-norm 
approximation. 


KRONECKER’S  THEOREM 


Let  {/in},  n  =  0, 1, . . be  a  causal  sequence  of  unit  impulse  responses.  Corresponding 
to  this  sequence,  we  associate  two  quantities:  the  2-transform,  or  symbol, 


H{z)  =  YJhnZ~n, 

n= 0 


and  the  infinite  Hankel  matrix  T h  =  1  <  i,£  <  oo,  or 


rw  = 


hi 

h.2 

hi 


hi 

hi 


(1) 


(2) 


Note  that  we  do  not  include  ho  in  the  definition  of  T#  in  (2)  simply  because  of  the  standard 
convention  in  systems  theory.  The  constant  ho  is  easy  to  determine  since  it  is  the  limit  of 
H(z)  as  2  tends  to  infinity. 

The  rank  of  the  matrix  Th  is  defined  by  the  number  of  linearly  independent  columns 
of  T;/;  or  equivalently,  it  is  the  dimension  of  the  range  of  T//  on  £2,  where  £2  denotes  the 
space  of  square- summable  sequences.  The  following  classical  result  due  to  Kronecker  can 
be  found  in  [4],  A  proof  is  included  in  the  Appendix  for  easy  reference. 


Theorem  1.  (Kronecker) 

The  infinite  matrix  T#  in  (2)  has  finite  rank  if  and  only  if  the  symbol 


OO 

'  (3) 

71=  1 


is  a  strictly  proper  rational  function  in  z;  that  is, 


Y,h*z 


P(z) 

Q(z) 


(4) 


where  P(z )  and  Q(z)  are  relatively  prime  polynomials  in  z  with  degree  ( P )  <  degree  ( Q ). 
Furthermore,  in  this  situation, 

rank (F h)  =  degree  ( Q ).  (5) 

In  fact,  if  T//  has  finite  rank  k ,  say,  then  the  first  k  columns  of  Th  are  linearly 
independent,  and  there  exist  k  constants  ciT. . .  ,c*,  such  that 
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7*+/  =  ^c,7;+/-i,  ^  =  1)2, ... . 


(6) 


i=i 


where  7,  denotes  the  ith  column  of  Th-  Furthermore,  the  finite  square  Hankel  matrix 


Hk 


hi  hi  ...  hk 
hi  .  hk-\ 


L  hk  .  h2k-i  J 


(7) 


which  is  a  finite  square  truncation  of  T/f ,  is  nonsingular.  For  more  details,  see  Lemma  2.1 
and  its  proof  in  [2]  or  the  Appendix  in  this  report. 


THE  ARMA  MODEL 

A  sequence  {/i„},  n  =  0, 1, . . .,  is  said  to  represent  an  ARMA  model  if  its  2-transform 
H(z)  defined  in  (1)  is  a  proper  rational  function  in  z  with  all  its  poles  lying  in  the  open 
unit  disc  \z\  <  1,  namely: 


H(z)  =  Y,h^~n  =  T 

n=0  1 


I  _  1TJ 


+  aiz  1  -| - +  aNZ~ 

bp zp  +  bi zp~1  -) - f-  bp  _  P(z ) 

zP  +  aizP-1  - 1- op  <5(r) 


(8) 


for  arbitrary  M,  N  >  0,  where  p  =  max(M,  Ar),  &Af+i  =••■  =  &„  =  0,  av+i  =  ■  ■  ■  =  ap  = 

0, 

degree(P)  <  degree(Q)  =  p,  (9) 

and 


H(z)  is  analytic  in  1*1  >  1-  (10) 

If  P(z)  and  Q(z )  have  no  common  factors,  we  say  that  the  representation  (8)  of  H(z)  is 
in  coprime  form.  If  (8)  is  in  coprime  form,  then  (10)  is  equivalent  to 

Q(z )  ^  0  for  all  \z\  >  1.  (11) 

We  remark  that  (10)  is  also  equivalent  to 

OO 

J>„l<°o-  (12) 

n=0 


In  the  following,  we  derive  an  algorithm  to  yield  the  feedback  sequence  {an}  and  feed¬ 
forward  sequence  {6n}  in  (8)  from  the  sequence  of  unit  impulse  responses  {An},  assuming 
that  its  2-transform  H(z)  is  known  to  be  an  ARMA  model. 
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Algorithm  I.  (ARMA  realization) 

(1°)  Solve  for  ai, . . . ,  ap  in  the  linear  matrix  equation: 


"  h\  h-2 

..  hp  i 

hp+i 

dp 

dp—\ 

+ 

hp- H 
hp+ 2 

'-^p . 

■  h2p- 1. 

.  ai  . 

-  h2p  - 

(13) 


°)  Compute  bo, bp  by  matrix-vector  multiplication  in: 

'1  0  .  O' 

ai  1  '•  : 

••  0 

dp  .........  ctj  1 

where  a\ , . . . ,  ap  have  been  computed  in  (13). 

Remarks:  From  (7),  with  k  =  p,  we  see  that  the  coefficient  matrix  of  the  matrix  equation 
(13)  (which  is  a  finite  Hankel  matrix)  is  nonsingular.  The  coefficient  matrix  in  (14)  (which 
is  a  finite  Toeplitz  matrix)  is  a  lower  triangular  matrix  with  unit  diagonal. 

Proof:  Multiplying  Q(z)  to  H(z)  in  (8),  we  have 

(ho  h\z  1  +  •  •  -)(zp  +  ai  zp  1  +  •  •  •  +  ap)  =  bozp  +  ■  ■  ■  +  bp.  (15) 

Hence,  equating  the  coefficients  of  zp,zp~1 ,z,  1,  we  have  (14);  and  equating  the  coef¬ 
ficients  of  z~\. . . ,  z~v ,  we  have  (13). 


ho 

hi 

L  hv\ 


(14) 


THE  HANKEL  NORM 


As  mentioned  in  the  Introduction,  we  must  assume  for  all  practical  purposes  that 
H(z)  is  not  a  rational  function.  Hence,  we  usually  start  with  a  sequence  {/i„}  of  unit 
impulse  responses  which  satisfies  ^  |h„|  <  oo.  In  fact,  in  most  situations,  we  only  have 
input/output  information  {un}/{vn}.  In  this  case,  we  write 


H(ej“) 


m _ 

]T)  ume~imw 


(16) 


and  numerically  compute  the  unit  impulse  responses  hn  by  applying  DFT  to  the  integral 


(17) 
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We  do  not  dwell  on  this,  since  it  is  expected  that  other  more  efficient  procedures  will 
be  introduced  in  our  forthcoming  report  [9] ,  but  only  remark  here  that  any  numerical 
approximation  method  for  determining  {hn}  from  (17)  produces  an  error,  which  again 
necessarily  results  in  a  non-rational  model  H(z ),  z  =  eJU\ 

Now’,  as  a  result  of  the  condition  ^  |/in |  <  oo,  we  have  a  bounded  function  H(z)  on 
the  unit  circle  |z|  =  1.  Hence,  we  may  write 


H(z)  =  Ha(z)  +  Ha(z), 


(18) 


where 


Ha(z) 


is  called  the  analytic  part  of  H(z)  and 

/  OO  \  OO 

H.(z)=  £  hnZ~n)  =Y.hnZ~n 


n=l 


(19) 


(20) 


is  called  the  singular  part  of  H(z).  Hence,  the  Hankel  matrix  T /;  corresponding  to  H(z). 
as  defined  in  (2),  is  uniquely  determined  by  the  singular  part  H,{z)  of  H(z),  while  the 
analytic  part  Ha(z)  of  H(z)  does  not  influence  T //  at  all.  Now,  under  the  hypothesis 
\hn  \  <  oo,  which  is  equivalent  to  BIBO  (bounded  input  bounded  output)  stability,  w'e 
may  conclude  that  T  //  is  a  bounded  linear  operator  on  the  space  £2  of  all  square- summable 
sequences,  with  the  operator  norm  defined  by 


l|r«||  =suP{||rHx||,2:  ||x||„=l},  (21) 

where,  for  x  =  (xi ,  x2,  ■ . .), 

11x11,2  =  (^T  |i, 

V.=! 

The  finiteness  of  ({T // 1|  follows  from  Nehari’s  theorem  to  be  stated  below;  and  it  will  be 
seen  that,  in  fact,  we  have 


oo  \  oo 

El'1"!2  <  IIChII  <  53 |A„|.  (23) 

n=l  /  n=l 

We  will  verify  (23)  after  we  state  Nehari’s  theorem;  but  at  this  point  it  should  be  noted 
that  ||r//||  is  the  spectral  radius  of  the  operator  Th,  or  equivalently,  the  largest  singular 
value  of  F  fj. 

The  operator  norm  of  Th  in  (21)  is  also  called  the  Hankel  norm  of  the  function  H(z  ) 
and  will  be  denoted  by 


l|ff(*)llr  =  l|rn||. 


Actually,  this  is  not  a  true  norm  but  only  a  semi-norm  since 


(24) 
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||tf(z)||r  =  ||tfU)-<7(*)llr  (25) 

for  any  function  g(z)  £  H°°.  Here  and  throughout,  H°°  denotes  the  space  of  all  bounded 
analytic  functions  in  |z|  <  1.  It  is  also  called  a  Handy  space.  The  importance  of  this 
observation  is  that  it  provides  a  motivation  of  the  celebrated  theorem  due  to  Nehari  [5]. 

Theorem  2.  (Nehari) 

Let  H(z)  be  a  bounded  function  defined  on  the  unit  circle  |z|  =  1.  Then 

PfO)l|r=  irf||/f(*)-s(r)||i-(W=„.  (20) 


We  now  derive  the  inequalities  in  (23).  First,  note  that  F n  —  T u,  •  By  (26),  with 
g(z )  =  0,  which  is  certainly  in  i7°°,  we  have 


l|r»li  =  \\Th.W  =  ||tf.(*)||r  <  ||H.(*)IU-(I,,«S1)  (27) 


oc 


oo 

<£  I'1"'- 

n  =  1 


This  verifies  the  upper  bound  in  (23).  To  verify  the  lower  bound,  let  e}  =  (1,0.0....). 
which  has  unit  ( 2  norm  so  that  from  the  definition  (21),  we  have 


llrwll  >  ||r#/ei H#a 

=  IK^i ,  -  -  •)!!/» 


(28) 


In  the  above  formulation,  we  have  also  established  that  the  Hankel  norm  lies  between  the 
L°°{\z\  =  1)  and  £2(|z|  =  1)  norms.  Recall  that 

||/lk°°(M=i)  =  sup  |/(s)|,  (29) 

l*l=i 

which  has  already  been  used  in  (26)  and  (27),  and  that 

ll/IU’iiii=i)  =  (2.  J*'  |/(e>»)|^  ’ .  (30) 

Indeed,  by  the  isometry  between  t2  and  L2(|z|  =  1)  we  have  the  lower  bound,  and  by  the 
application  of  Nehari's  theorem  as  in  (27)  we  have  the  upper  bound;  that  is, 

l|tf*IU*(|*|=i)  -  IWIr  <  ||^IU«(P|=1)-  (31) 
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SINGULAR  VALUES  AND  SINGULAR  VECTORS 


Let  T h  be  an  infinite  Hankel  matrix  as  defined  in  (2).  Under  the  assumption  that 
V  |/in|  <  oo,  we  have  seen  that  Th  can  be  considered  as  a  bounded  linear  operator  on  t.2 . 
Fbr  simplicity,  we  will  call  T //  a  Hankel  operator.  Since  all  Hankel  operators  are  symmetric, 
the  adjoint  of  T h  is  its  complex  conjugate 


Th  = 


h 


/13 


(32) 


In  particular,  for  a  sequence  of  real- valued  unit  impulse  responses,  Th  is  self-adjoint.  In 

any  case,  r  h  F 11  is  non-negative  definite  and  has  non- negative  eigenvalues  s2 ,  s\ . with 

sj  >  ^2  >  •  •  •  >  0,  which  are  arranged  in  nonincreasing  order,  listing  ah  multiple  occur¬ 
rences.  Let  jr // 1  denote  the  non-negative  square-root  of  ThTh-  Note  that  although  |T // 1 
is  Hermitian,  it  may  not  be  Hankel;  but  this  is  not  important  in  the  following  discussions. 
Clearly,  the  eigenvalues  of  jT// 1  are  sj,S2,...  .  Corresponding  to  each  s„,  let  x„  be  an 
eigenvector  of  [T// 1;  that  is,  xn  ^  0,  x„  G  £2,  and 

|rw|x„  =  s„xn.  (33) 

Let  U  be  an  unitary  operator  on  £2  such  that 


(cf.  Ref.  10).  Then  by  defining 


r«  =  tf|r«| 


Vrx  =  UX 


we  see  that 


Indeed,  it  is  clear  that 


F/i^n  — 

F  H^]n  =  Sn£n. 


r//(n  =  u\rH\U 


—  SnF  —  Sn^ln 


(34) 


(35) 


(36) 


(37) 


and 


rHvn  =  r„utn 

=  (uTrH)Ttn 

=  \tn\tu  =  |rw|xn 


(3S) 


Chui 


We  will  call  s„  a  singular  va'ue  (or  5-number)  of  Th  and  will  call  the  pair  (^n^n)  a 
singular  vector  pair  (or  Schmidt  pair)  of  Th  relative  to  s„.  Let  us  assume,  without  loss 
of  generality,  that  the  eigenvectors  {xi,x2,...}  of  |T h |  corresponding  to  the  eigenvalues 
{si,52,  . . .}  are  orthonormalized,  so  that  the  infinite  matrix 

V  =  [xix2  . . .],  (39) 

with  xn  as  its  nth  column,  is  a  unitary  operator  in  t2.  Then,  by  the  definition  of  eigenvalues 
and  eigenvectors,  we  have 


|rw|v  =  vs, 

where 


(40) 


S  =  diag(si,s2,...,)  = 


•Si  O 

S2 


O 


(41) 


is  the  diagonal  matrix  of  the  singular  values  of  Th-  Hence,  it  follows  from  (34)  that 


rH  =  (UV)  SV*  (42) 

where  V*  is  the  adjoint  (or  complex  conjugate  of  the  transpose)  of  V,  and  both  (UV)  and 
V*  axe  unitary  operators.  The  decomposition  (42)  is  called  a  sin(  r  value  decomposition 
of  Th-  (For  finite  matrices,  see  Ref.  11.1  Another  formulation  ol  (42)  is  the  so-called 
Schmidt  series  representation 


OO 

Fh  =  ]Psn?7n{*,  (43) 

n=l 

which  follows  from  (42)  and  (35);  where  again,  denotes  the  complex  conjugate  of  the 
transpose  of  the  vector  so  that  rjn£*  is  a  rank-1  operator  on  £2.  That  is,  Th  is  the 
(strong)  limit  of  the  sequence  of  finite  sums 

p 

(44) 

n=l 


(which  are  rank  p  operators),  as  p  — >  oo.  “Strongness”  here  means  convergence  in  the 
operator  norm  defined  in  (26).  However,  although  the  operators  in  (44)  are  finite-rank 
operators  that  approximate  Th,  they  are  not  Hankel  operators,  ana  consequently,  are 
useless  for  our  purpose  of  finding  rational  models  for  H(z )  (see  Theorem  1). 


AAK  THEORY 

Let  H(z )  and  T //  be  defined  as  in  (1)  and  (2)  where  the  sequer  e  {/in}  of  unit  impulse 
responses  is  assumed  to  be  in  U ;  that  is, 

OO 

^2  !M  <  oo.  (45) 

n=l 
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Let  si  >  52  >  •  •  •  >  0  denote  the  singular  values  of  T//  with  corresponding  singular 
vector  pairs  (6, r?i), (£2, 1)2)1  ■  •  •>  respectively.  Since  Si  is  the  largest  singular  value,  it  is 
the  spectral  radius  of  T,  so  that  we  have  sj  =  ||r 7/ 1|.  Hence,  as  a  consequence  of  Nehari’s 
theorem  (cf.  Theorem  2),  we  have 

si  =  IM  =  || H(z)  -  5(2)||i=o(|z|=1).  (46) 

That  is,  the  first  singular  value  si  yields  the  exact  measurement  of  the  error  of  uniform 
approximation  on  |z|  =  1  of  H(z)  from  the  Hardy  space  H°°.  That  is  why  Nehari’s  result 
is  considered  to  be  the  starting  point  of  the  area  of  ff°°-control  (cf.  Ref.  12).  To  yield 
rational  models,  we  introduce  the  notation: 

ftp  =  {g(fj:  de&P(z)  <  deiQ(z)  <  Pi  Q(z)  ±  0  for  lz\  >  1  j  (47) 

where  P,  Q  are  polynomials.  That  is,  TZp  is  the  collection  of  all  stable  proper  rational 
functions  of  degree  at  most  p.  Also,  let 

Hoo  =Hoo  +  1lS)  (48) 

denote  the  collection  of  all  functions  f(z)  =  g(z)  -+■  r(z)  where  g(z)  6  H°°  and  r(z)  E  K,p. 
To  unify  notations,  we  set 

H°°  =  P0°°,  TVq  =  {0}.  (49) 

To  pose  the  approximation  problem  in  terms  of  operators,  let  Gp  denote  the  collection  of 
all  bounded  Hankel  operators  with  rank  <  p.  Hence,  Go  consists  only  of  the  zero  matrix. 
By  Kronecker’s  theorem  (cf.  Theorem  1)  and  the  equivalence  of  (11)  to  the  boundedness 
of  the  corresponding  Hankel  operators,  we  may  identify  7lp  with  Gp. 

Hence,  Nehari’s  theorem  (cf.  Theorem  2)  formulated  as  in  (46)  can  be  written  as: 

Sl  =  A€Go  ^  =  ~  5(Z)Hf'00(ld=l)*  (5°) 

Here,  since  Go  =  {0},  the  first  equality  only  says  Si  =  ||r h\\-  The  celebrated  result  of 
AAK  in  [3]  is  to  generalize  (50)  from  p  =  0  to  an  arbitrary  positive  integer  p. 

Theorem  3.  (AAK) 

Let  H(z)  be  any  function  in  L°°(|z|  =  1)  and  p  be  any  non-negative  integer.  Then 

VH  =  j**  llr*  ~  All  =  j|^00  I \H{z)  -  0(z)||L~(M=1).  (51) 


Hence,  the  (p  +  l)st  singular  value  of  Fn  gives  the  exact  measurement  of  the  distance 
in  uniform  norm  on  |z|  =  1  of  H(z)  from  H£°.  In  addition,  since  Th  is  a  compact  operator 
under  the  assumption  (45),  we  have 

sp+i  — ♦  0  as  p  — *  00  (52) 

(cf.  Ref.  13).  That  is,  we  can  approximate  H(z)  as  close  as  we  wish  from  H£°.  So,  if 
the  Hankel  norm  defined  in  (24)  is  used,  then  in  view  of  (25),  the  P°°  or  analytic  part  of 


9 


Chui 


the  approximant  is  irrelevant.  Thus,  by  taking  the  singular  part  [  ]a,  we  obtain  a  stable 
rational  model  of  H(z) 

To  be  more  specific,  we  exhibit  the  formulation  of  the  best  Hankel-norm  approximant. 
The  following  notations  are  needed.  Let 


fp+i  =  («i 


(p+i)  „>+ 1) 


u 


,...),  T)p+1=(v\ 


(p+1)  ,.(P+1) 


(53) 


be  a  singular  vector  pair  of  Fh  corresponding  to  the  ( p  +  l)st  singular  value  sp+i  of  rn- 
Define  the  analytic  and  singular  functions: 


»  (54) 

.  i=  1 

respectively,  and  set 

9P(Z)  =  H(z)  -  sp+ i^TT-  (55) 

In  the  AAK  paper  [3],  it  is  proved  that  gp(z )  is  in  H£°  and  solves  the  optimal  approximation 
problem  (51)  in  the  sense  that 


sp+ i  =  llrH  “  r*J  =  \\H(z)  -  9p{z)\\l~(\z\=\).  (56) 

As  we  remarked  above,  since 


T9p  =rlM- 

where  [<)p]s  denotes  the  singular  part  of  gp ,  we  may  conclude  that 


(57) 


sP+i  =  || H  -  [gp]s||r  =  inf  || H{z)  -  g(z)||r;  (58) 

9 

and  this  means  that  the  singular  part  of  gp{z),  defined  by 

&p(z)  =  [ffp(2)]*»  (59) 

provides  an  optimal  stable  rational  approximation  of  H (z)  from  77* .  The  proof  of  this  result 
is  very  deep;  but  if  H(z)  itself  is  also  a  rational  model,  a  simpler  derivation  is  accessible. 
We  will  outline  this  simpler  version  in  the  next  section.  The  reason  for  doing  so  is  that 
S.Y.  Kung’s  model  reduction  algorithm,  which  will  be  discussed  in  this  report,  is  based  on 
this  derivation. 


DISCUSSION  OF  AAK’S  THEOREM 

We  will  outline  the  proof  of  AAK’s  theorem  (cf.  Theorem  3)  for  the  special  case 
where  H(z)  =  ^  hnz~n  has  real  coefficients  and  the  singular  part  [J7(z)]j  =  hnz~n 

of  H(z)  is  a  strictly  proper  rational  function,  namely: 


10 


NRL  Memorandum  Report  6689 


[H(z)]a  =  Ha(z)  = 


Pm-i(z) 

Qm(z) 


with  degree  (Pm- i)  <  degree  (Qm)  —  M,  and  Pm~i(z),Qm(z)  being  relatively  prime. 
Note  that  Ha(z)  6  TVM  since  T#  is  a  bounded  Hankel  operator. 

Let  be  the  (necessarily  non-zero)  eigenvalues  of  the  read  rank-Af  Hankel  matrix 

T ff.  We  arrange  these  eigenvalues  in  such  a  way  that  |A,j  =  s,;  hence, 

|Ai|  >  |A2|  >  •  ■  >  |Am|  >  0,  (61) 

where  Si  is  the  ith  singular  value  of  Th.  (Note  that  since  T#  is  real,  ThTh  =  ;  and 

this  gives  s?  =  A^.)  Let  x*  be  the  eigenvector  of  Fh  corresponding  to  Aj,  i  =  1 
Then  since  s,  =  |A,|  =  A,(sgn  A,),  we  have 


{Th6  =  siVt 
Thtu  =  S&  ' 


where  £,  =  x^  and  rji  =  (sgn  A,)xj.  That  is,  (x^,(sgn  A ,)x;)  =  is  a  singular  vector 

pair  of  F h  relative  to  the  singular  value  s*  =  |Aj|. 

Let  p  <  M,  and  our  goal  be  to  replace  Ha(z)  =  Pm-i(z)/Qm(z)  by  a  rational  function 
in  71^.  The  notations  introduced  in  (53)  and  (54)  axe  now  used. 

Step  1.  We  claim  that 

H(z)  =  Ha(z)£++i(z)  -  sp+1T]~+i(z)  (63) 

is  analytic  in  |z|  <  1;  that  is  [ff(z)]a  =  0.  To  see  this,  first  observe  that 

sp+iVp+i  ~  Ap+i^p+i  (64) 


so  that 


Sp+iVp+ii*)  =  Ap+i  ^u^+1)z  \ 


Second,  it  is  clear  that 


(p+l)  /— t— 1 
/  2 


Y  Y  h'u(t 

«=i  t= i 


= E  E  + £  E  A,-i+1^+y  ,2- 


i=]  t=i 


i=i  i=i 


where  ht  =  0  for  t  <  0.  Finally,  from  Th(p+i  =  Ap+1£p+i,  we  have 

OO 

£fc.+,_1u<'+’,  =  A,,+1u<''+,|,  «  =  1,2, ...  . 
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Hence,  we  may  conclude,  by  applying  (64),  (65),  (66),  and  (67)  consecutively,  that 


(68) 


=  !ff.wC+1w  -  ap+,  Y, 


:=i 


EE  EE  hi—  i+lu/-fi  2  ^  ^  -^p-Huj 


(?+l)z-t 


Li=l  t=  l 

oo  oo 


1=1  t=\ 


t=l 


YlYlh‘-'+iu 


(p+l)J-l 
/+1  z 


L«=i  (=1 

OO  OO 

e  e  Mr1* 

Li=i  <’=i+i 


(p+1)  — i  — 1 


=  0. 


Step  2.  We  claim  that  the  function 

B(z)  =  Qm(2)Vp+A2)  (69) 

is  a  polynomial  with  real  coefficients  and  of  degree  <  M  —  1. 

It  is  clear  that  B(z )  has  real  coefficients.  We  first  observe  that  since  Pm-\(z)  and 
£)j+i(z)  have  only  positive  powers,  we  have 


[pm-iM«;+i(*)].  =  °- 

Hence,  it  follows  from  Step  1  (cf.  Eq.  (63))  that 

[B(z)\s  =  [QM(z)T)~+i(z)}s 

■  [Pm-i(z)Zp+i (2)  -  sp+iQm (z)rip+i («)]4 


(70) 


(71) 


sp+i 

1 

■Sp+l 

1 

•Sp+l 


[Qm{z){Hs{z)Z++1{z)  -  Vhi?p+1(z)}L 
[Qm(z)H(z)\,  =  0. 


That  is,  B{z)  is  analytic  and  by  the  definition  (69),  it  must  be  a  polynomial.  Now,  since 


lim 


B(z) 


]  =  W*) =  °> 


PI— OO  Qm(z) 

we  have  degree  (B)  <  degree  (Qm),  so  that  the  degree  of  B(z)  is  at  most  M  —  1. 


(72) 
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Step  3.  Let  B*(z)  and  Qm(z)  be  the  reciprocal  polynomials  relative  to  B(z )  and  Qm(z), 
respectively;  that  is, 

B*(z)  =  zm~1b(±)  and  Q*M(z)  =  zMQN  ( ±) .  (73) 


We  claim  that  the  function 


C(z)  = 


Pm-\{z)B*(z)  -  Xp+1Q*M(z)B(z) 


w  Qm{z) 

is  a  polynomial  with  real  coefficients  and  of  degree  <  M  —  1. 

Indeed,  since  Pm-\{z)B*(z )  —  A p+iQ^f(z)B(z)  is  a  polynomial  of  degree  <  2 M  —  1, 
we  may  use  partial  fractions  to  write 

PU-,(z)B‘(z)-\+IQ\,(t)B(z)  C(z)  D(z) 

Qm(z)QU 2)  C«W  '  1 

for  some  polynomials  C(z)  and  D(z)  of  degree  <  M  —  1.  However,  since  Q m{z)  has  all  its 
zeros  in  \z\  <  1,  its  reciprocal  polynomial  Q*m(z)  has  all  its  zeros  in  \z\  >  1,  so  that 


Qm(z) 


Qm(\ ) 


Y,*'2 1 


for  some  <?o,  q\ ,  •  •  •  .  Hence,  if  degree  ( D )  <  M  —  1,  we  have 


IQm(z)1*  [  ^ 


fD(z)Y 


D{z) 

q'Z  Quiz)' 


On  the  other  hand,  since  Q*M(z)  is  analytic  in  \z\  <  1,  we  have 

[£W1  -o  (78 

Now,  in  view  of  Step  1,  namely:  \H(z)\a  =  0,  we  may  apply  (69)  in  Step  2  to  verify  that 

Pm-i{z)B*(z)  —  \p+iQ*m(z)B(z)'\  _ 

Qm{z)Q'm{z)  J,  U’  ^  ■ 

so  that  by  combining  (77),  (78),  and  (79)  in  (75),  we  have  D(z )  =  0.  This  verifies  (74). 
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Step  4.  Let  Rp(z)  be  defined  as  in  (59).  We  claim  that 

-  ,  C(z)  ' 

Rp(z)~  [b-(2)].‘ 


Indeed,  from  (74),  we  have 


C(z )  _  rrf  ,  ,  B(z)/Qm(z) 

B*(z)~H{z)  Xp+1  B*(z)/Q*m(zY 


Here,  from  the  relations  (69)  and  (64),  it  can  be  shown  that 

'  B(z)/Qm(z )  Vp+iif)  (82) 

^’b-m/qw*)],  ”+1  L^+.Wj,' 

Therefore,  in  view  of  the  definitions  in  (55)  and  (59),  we  have  obtained  (80). 

Step  5.  We  can  now  complete  the  proof  of  the  AAK  theorem  (cf.  Theorem  3),  for  the 
special  case  of  real  and  finite-rank  Hankel  operators  T//,  by  verifying  (56),  or  equivalently: 

Sp+i  =  ||IVr  -r£pll-  (83) 

First,  we  remark  that  VP+i{z) / £p+i(z)  's  a  constant  c  multiple  of  a  Blaschke  product  with 
|c|  =  1,  so  that 

=  1.  (84) 

£p+ 1  (  2  )  £,<x>(|z|  =  l) 

Therefore,  from  (55)  and  (31),  we  have 

l|r»-rs||  =  ||rH-r,,||  m 

=  \\H  -$Jr  =  V-i  I??1  r 

Vh  1 

-Sp+1  7+17\  ~  p+ l' 

Cp+n2;  o*i=i) 

On  the  other  hand,  since 


[(^»(z)  "  ^p(5))fp+i(2)]»  ~  t+  /7\ 


=  »,+.  fe^kf+.wl  =  vulval- 

tp+Az)  J, 

=  ■Sp  +  l  T)p4-\ (z)i 
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we  have 

rH-fi/p+i  ^p+i^7p+i»  (87) 

so  that 

■*■  u_  ^*7p+i  ^p+i(p+i-  (88) 

That  is,  (.sp+1,(£p+i,77p+1))  is  also  a  (singular  value,  singular  vectors)  pair  of  j  . 
Hence,  sp-f i  does  not  exceed  its  spectral  radius,  namely: 

*,«<lir„_Sfll  =  ||r«-rgj|.  (89) 

Therefore,  by  combining  (85)  and  (89),  we  have  proved  (83). 

Remark.  Prom  (80),  it  seems  that  Rp(z )  is  a  rational  function  in  7ZSM_1.  Actually,  we 

have  now  proved  that  Rp(z )  is  in  where  p  <  M  —  1.  This  means  that  only  the  p  zeros 

of  B*(z)  that  lie  in  \z\  <  1  are  used  to  yield  Rp.  The  results  in  Steps  2,  3,  and  4  are  useful 
in  the  following  discussion  of  Kung’s  model  reduction  algorithm. 


KUNG’S  MODEL  REDUCTION  ALGORITHM 

Let  H„(z )  =  Pm-i(z)/Qm(z)  be  defined  as  in  (60)  where  Pm-i(z)  and  Qm{z)  are 
coprime  polynomials  with  degree  (Pm- i)  <  degree  (Qm)  =  M.  Let  Q*M{z)  be  the  recip¬ 
rocal  polynomial  relative  to  QM(z)  as  defined  in  (73).  Since  Hs(z)  is  in  7lsM,  all  the  M 
zeros  of  Qm(z)  lie  in  \z\  <  1,  so  that  all  the  M  zeros  of  Q*M(z)  lie  in  |z|  >  1.  That  is, 

^hl  =  ,0  +  „z-1+,2z-2  +  ...  (90) 

is  the  reciprocal  of  a  finite  Blaschke  product. 

We  need  the  following  notations: 


Hm  = 


Km  = 


hi  h2 
/l2 


hm 

hM+i 


Ha  = 


Km  .  h2M-\  J 

9m  9m-i  •  •  •  9i 
9M+1  9M  92 

92M-1  92M-2  •••  9M 

0  .  0 


hi 


0  hi  ...  hM-i 


(91) 


(92) 


(93) 
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r  9o  o  ...  0-1 


K 


A  = 


9i 


L9M-1 


9i  9o  -I 


(94) 


where  Hm,  H&.  axe  finite  Hankel  matrices  and  Km,  K&  axe  finite  Toeplitz  matrices.  (Note 
that  the  notations  in  (90)  and  (7)  agree.) 


Algorithm  II.  (Rung) 

(1 i)  Solve  the  generalized  eigenvalue  problem 

{Hm  —  \Km)c\  -  0  (95) 

for  {A,,  q'},  where  {q1 , . . . ,  qM}  are  linearly  independent  and 

|Ai|>|A2|>...>|Am|.  (96) 


(2°)  For  1  <  p  <  M,  perform  the  matrix-vector  multiplication: 


rp+1  =(HA-\p+1KA)qr+'. 


(3°)  Set 


and  compute  the  singular  part 


rp+1  =  (b0,. . .  ,bM~\) 
qp+1  =  (a0,...,QM-i)’ 


hozM  1  +  •  •  •  +  5m— i 
ao  +  aiz  +  •  •  ■  +  clm -\zM~ 1 


(97) 


(98) 


(99) 


(In  this  computation,  factorize  ao  +  oiz  -f-  •  •  •  +  clm-\zm~ 1  into  M  —  1  linear  factors  and 
aiscard  those  terms  whose  roots  do  not  He  in  \z\  <  1.  We  note  that  exactly  p  linear  factors, 
counting  multiplicities,  axe  retained.  Now  use  partial  fraction  expansion  to  determine 
Rp(z).) 

Proof.  Let 


’  hM+ 1 

hM+2  ••• 

h2M 

hM+2 

. 

h2M+l 

,  (100) 

-  h2M 

. 

hiM-i  - 

<l2M 

92A/-1 

9m+ 1  ■ 

92M  +  1 

92M 

9M+  2 

(101) 

-  93  At  —  1 

93M-2  ••• 

92  Af  - 
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where  H M  is  a  finite  Hankel  matrix  and  KM  a  finite  Toeplitz  matrix.  Set 

C(z)  =  bozM  1  +  •  •  •  -f  b\i- 1 
B(z)  =  ao  1  4-  ■  •  ■  +  OLM- 1 

so  that 

B  (2)  =  ao  +  ai z  -f-  •  •  •  -(-  1 

Then  it  is  straightforward  to  check  that  the  identity  (74)  in  Step  3  in  the  proof  of  AAK’s 
theorem  is  equivalent  to  the  matrix  system 


(102) 


(103) 


H  a 

Hm 


a  0 


LajW-i  J 


-  A 


p+ 1 


Ka 

Km 

km 


Oq 


Lqm-i  J 


&AJ-1 

0 

0 

0 


Also,  (104)  is  equivalent  to 

'[i/A-Ap+1KA]q^+1)=r(P+I) 
-  [tfM-A,+1tfM]qO>+1>  =  0 
[ffM  -  Ap+1ATM]q(P+1)  =  0 


(104) 


(105) 


where  the  notations  in  (98)  are  used.  The  second  equation  in  (105)  is  Step  (1°)  in  Algo¬ 
rithm  II,  while  the  first  equation  in  (105)  is  Step  (2°)  in  the  algorithm.  Note  that  the  third 
equation  in  (105)  is  a  consequence  of  the  second  equation  by  the  Caley-Hamilton  theorem 

in  matrix  theory.  Now,  the  result  Rp(z)  computed  in  Step  (3°)  of  the  algorithm  agrees 
with  Rp(z)  in  (80);  and  hence,  by  AAK’s  theorem  (cf.  argument  in  Step  5  in  the  above 
section),  we  may  conclude  that  Rp{z )  satisfies  (83). 


REFORMULATION  OF  OPTIMAL  HANKEL-NORM  RATIONAL  MODEL 

We  have  seen  that  Toeplitz  matrices  also  occur  in  the  computation  of  the  rational 
model  Rp(z )  in  Kung’s  algorithm.  In  fact,  a  Toeplitz  matrix  associated  with  the  transfer 
function  H(z )  can  be  used  to  describe  Rp(z).  Let 


17 


Chui 


'0  h\  /12  h$ 

0  0  hi 

0  0  0  hi 


(106) 


be  an  upper  triangular  infinite  Toeplitz  matrix  with  zero  main  diagonal.  As  in  (61),  let 
{Afjjlj  denote  the  (necessarily  non-zero)  eigenvalues  of  the  rank-M  Hankel  operator  T h 
such  that  | Ai  j  =  s,  is  the  ith  singular  value  of  T//,  t  =  1, . . .  ,  M.  As  before,  let 

m  =  =  (107) 

denote  an  eigenvector  of  T  h  corresponding  to  the  eigenvalue  A,,  so  that  (£,,171),  with 
T],  =  (sgn  A,)xj,  is  a  singular  vector  pair  of  T h  corresponding  to  the  singular  value  st  =  |A,| 

(cf.  Eq.  (62)).  Also,  let  gp(z)  be  as  in  (55)  and  Rp(z)  =  [gP(z)]a  be  the  optimal  solution 
in  best  Hankel-norm  approximation  of  H(z)  from  1Z*.  Then  we  have  the  following  result. 

Theorem  4.  Let  H(z )  be  a  transfer  function  with  real  coefficients  such  that  its  singular 
part  Ha(z)  is  given  as  in  (60).  Then 


RP(z)  = 


(1,  2,  Z  •  ■  .)  TjfXp.|-i 

(l,z,z2,...)  -Xp+1 


(108) 


Proof.  From  the  equation  T,tfXp+1  =  Ap+1xp+],  we  have 

OO 

Eh-  .  —  \  D  _  I  O 

"■i+t— 1 u ,  Ap-)_i  Uf  ,  {  —  1,2,..., 

1=1 

and  hence,  by  (65),  it  follows  that 


(109) 


OO  /  OO 


sP+iVP+i(z)  =  Vh 


b- 


(110) 


This  yields: 


Hs{z)tp+i(z)  -sp+iT7p+i(z) 


(in) 


OO  /  OO 


-51  ( ^hi+t~iut' 

i=  1  \/=l 


(pH-1) 


00/00 


=E(  £ 

i=  1  \/=t+l  / 

=(l,z,z2,...)  T/xp+1. 

Since  £p+1(z)  =  (1,  z,  z2, . . .)  •  xp+1 ,  we  have 
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rj  r  \  _  Vp+li2)  _  (E  Z,  Z2,  ■  •  •)  '  TfXp  +  ] 

ns{Z)  Sp+1  ,  ■ ,  2  \  v  ’ 

Cp+iv-2)  yi-i2!2  xp+i 

so  that  by  (59)  and  (55),  which  was  proved  by  AAK  [3]  in  general  rind  in  this  report  for 
our  special  case,  we  have  obtained  (108). 


TRUNCATED  HANKEL  OPERATORS 


To  apply  Rung’s  algorithm,  it  is  essential  to  start  with  an  ARMA  model.  However,  as 
we  mentioned  before,  in  all  practiced  purposes  in  underwater  acoustic  signal  processing,  the 
sequence  {hn}  of  unit  impulse  responses  which  may  be  computed  from  input /output  mea¬ 
surements  (cf.  Eqs.  (16)  and  (T7),  for  instance),  is  both  physically  and  numerically  noisy. 
Hence,  the  measured  transfer  lunction  H(z)  cannot  be  ARMA  or  the  Hankel  operator  Y h 
has  infinite  rank.  In  a  recent  work  [7],  the  truncated  (infinite)  Hankel  operators 


h\  h 2  ...  hn  0 

/?2 

hn 

.  o  O. 


(112) 


are  introduced  and  optimal  Hankel-norm  approximants  R"(z)  of 

n 

Hn(z)=  k'Z~'  (H3) 

i=  —  oo 


are  used  to  replace  the  optimal  Hankel-norm  approximant  Rv{z)  of 

OC 

H(z)=  Y.  h 

*=— OO 


(114) 


Recall  that  this  means: 

'  *i+.  =  -  A"WUr  =  „,iyf  II^.W  -  A(*-)llr 

\\r 

'  *P+1  =  ll»M  -  Wllr  =  „  inf  Jff(z)  -  K(*)||r 

where  s(jn)  >  >  •  ■  •  are  the  singular  values  of  Th„  =  T It  is  at  least  intuitively  clear 

that 


SP+ 1  — *  -Sp+i 


as  n  — >  oo; 


(116) 
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and  in  [8],  we  have  proved  much  more. 

Theorem  5.  For  each  p ,  there  exists  a  positive  constant  Cp,  such  that 


£  I'**! 


1/2 


+  (n  +  1)  (  ^  1^*1 

\k=n+ 1 


oc  /  oc 


+  Y 

i  =  n  +  2  \j  =  k 


(117) 


Hence,  if  { /*„ }  decays  to  zero  very  rapidly  (as  is  the  case  in  most  underwater  acoustic 
experiments),  we  know  that  Rp(z)  converges  uniformly  on  |r|  =  1  very  rapidly  to  the 
optimal  Hankel-norm  rational  model  Rp(z).  In  the  next  section,  we  will  give  an  efficient 
algorithm  for  computing  R”(z). 


ALGORITHM  FOR  OPTIMAL  RATIONAL  MODEL 

Let  {//„},  n  =  1,2,. . .,  be  a  sequence  of  unit  impulse  responses  that  satisfies 

OC- 

Y'.  \hn  |  <  oo.  (118) 

n=  1 


(Note  that  H(z)  =  ^  h„z  "  is  not  necessarily  a  rational  model.)  Given  a  tolerance  f  >  0. 
we  now  give  an  algorithm  for  computing  a  stable  (strictly  proper)  rational  function  Rp(z) 
such  that 


\\H(z)-Rp(z)\\r  <e.  (119) 

The  degree  p  will  be  chosen  to  be  the  smallest  possible,  under  the  limitation  of  this  method. 

Algorithm  III. 

(1°)  Choose  an  £j,  0  <  £\  <  £  ( e.g .  £\  =  ),  and  the  smallest  positive  integer  A/j, 

such  that 

OO 

Y  \hn\<£u  (120) 

n=Af,  +  l 

where  h\i}  /  0. 

(2°)  Find  the  singular  values 


(i) 


>  s2 


(i) 


>  >  -«M, 


>  0 


(121) 
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of  the  finite  Hankel  matrix 


■p(  i) _ 

1  H  ~ 


h  i  h2  •  •  •  h\]  i 

h“i 


and  choose  the  smallest  positive  integer  p\  <  M\  such  that 


4,  <e-£i. 


(123) 


( One  can  find  the  eigenvalues  of  and  then  take  their  absolute  values  to  determine  the 
singular  values.) 

(3°)  If  p i  does  not  exist  (i.e.,  if  4V,  >£  —  £])•>  return  to  (1°)  with  £j  replaced  bv 


£i 

£2  =  T 


and  Mj  replaced  by  M2,  so  chosen  that 


^2  \hn\<£-2  = 


n—  Mz  +  l 


and  ^  0.  M2  >  M\. 

Repeat  (2°)  by  finding  the  singular  values 


42>  >42>  >0 


of  the  Hankel  matrix 


p(2)  _ 

1  H  ~ 


hi  h2  ...  h\i7 

h2 


and  choosing  the  smallest  p2  <  M2  such  that 


4v  <  £  -  ^ 


(128) 


(4°)  If  this  fails  (i.e.  if  s ^  >  £  —  £2),  repeat  (3°)  with 


,  __  £2 
3  ~  T  ’ 


(129) 
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etc.  Suppose  that  the  first  time  this  procedure  succeeds  is  at  the  kth  iteration.  That  is,  k 
is  the  smallest  integer  such  that 

s(m\  <  e  -  £k  =  £  ~  1 


where 


Ak) 


s(k) 

sMk 


2k~1  £l’ 

(130) 

>  0 

(131) 

are  the  singular  values  of  the  M k-dimensional  square  Hankel  matrix 


r(*)  - 
1  H  ~ 


[  hi  h2  ...  hMk 
h2 


h\ik 


o 


with  h\ik  /  0,  and 


o°  1 

Y  Ife«l  <**  =  2*^*1  • 


Let  pk  he  the  smallest  integer  such  that  <  £  —  £k,  and  set 

P  =  Pk  - 1; 


so  that 


0  <  —  S<*>  <  E 

U  <  sMk  -  5p+l  —  SPk  -  £ 


2^£i' 


(5°)  Compute  the  eigenvalue-eigenvector  pair 

(Vt-i,q’’+1) 

ofr';1,  where  IVnl  =  s‘+V 

(6°)  Perform  the  matrix-vector  multiplication 

0  . 

I 

,j>+i 


0 

hi 


L0  h\  ...  h\ik-i 


(132) 


(133) 


(134) 

(135) 


(136) 


(137) 


(7°)  Write: 


{ 


|  rp+1  =  (b0, . . . ,  b\ik-\ ) 

lp+1  =  (ao, 


(138) 
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and  compute  the  singular  part 


bpzMk  + 

do  +  di  Z  -| - 1-  OLMk- \ZMk~1 


(139) 


by  using  partial  fractions,  retaining  only  the  p  poles  in  |z|  <  1,  counting  multipUcities. 
Then  Rp(z )  (E  7Z*  and  satisfies 

\\H(z)~  Rp(z)\\r  <£.  (140) 


Proof.  We  first  remark  that  det  =  (  —  l)Mkh^kk  ^  0  so  that  all  the  singular  values 

•s,  ,i  =  l,...  ,Mk,  are  positive  (nonzero).  From  a  well-known  result  in  operator  theory 
(cf.  [13]),  we  have 


,  max  |5, -s[fc)|  <  ||r// -r^H.  (141) 

1  <t<Mk 

On  the  other  hand,  it  follows  from  the  definition  of  rff  and  the  second  inequality  in  (23) 
that 

OO 

iir«  -  r^>||  <  £  im- 

n=A/*  +  i 

Hence,  by  combining  (141),  (142),  and  (133),  we  obtain 

i  (*■)!  ^  1 

max  \s,  -  s,  <  — r — ~£\  ■ 

1<,<M„  1  11  2*-1 

Now,  since  {hn}  satisfies  (118),  we  have  sn  — ►  0  by  (52);  so  that  in  view  of  (143)  (where  Mk 
necessarily  tends  to  infinity,  as  k  — >  oo),  it  follows  that  (130)  is  satisfied  for  all  sufficiently 
large  values  of  k.  This  shows  that  the  iteration  procedure  (l°)-(4°)  converges. 

Next,  from  the  definition  of  in  (132),  we  note  that  the  corresponding  rational 
symbol  of  is 


(142) 

(143) 


Mk 

tf*(z)  =  ^/il2-‘ 

i=l 


hiZAfk  1  +•••-(-  /lM, 

jMk 


P(z) 

<?(*)’ 


j 

(144) 


where  degree  (P)  <  degree  ( Q )  =  Mk  and  Q{z )  =  zMk.  Since  the  reciprocal  polynomial 
of  Q(z)  is  Q*(z )  =  1,  we  have 


Q*(z)  -Mk 

Q(z) 

Hence,  using  the  notations  from  (92)  and  (94)  with  M  replaced  by  Mk,  we  have: 


(145) 


KMk  = 


1 

o 


o 

1 


,*A  =  [0] 


(146) 
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Also,  again  from  the  definition  of  in  (132),  using  the  notations  from  (91)  and  (93) 
with  M  —  Mk,  we  have: 


and 


This  yields 


Ha  = 


Hm  =  r(, *> 

)  .  0 

hi 

L  0  h  i  ...  h  \fk  _  ] 


f  (H\ik  -  \KMk) q  =  r(Hfc)q  -  Aq 
\  (Ha  -  AA'A)q  =  Ha<1. 

Hence,  Step  (1°)  in  Kung’s  algorithm  (i.e.  Algorithm  II)  is  equivalent  to 

r(*)nP+1  —  A  . ,np+1 

and  Step  (2°)  in  Kung’s  algorithm  (i.e.  Algorithm  II)  is  equivalent  to 


(147) 


(148) 


(149) 


(150) 


rp+1  =  tfAqp+1.  (151) 

Since  (150)  is  the  same  as  Step  (5°)  and  (151)  is  the  same  as  Step  (6°)  in  Algorithm  III. 
we  have  proved  that  Rp(z),  as  defined  in  (139),  satisfies 

||H,(*)-Jfc(*)||r=#,„  (152) 

where  Hk(z)  is  defined  in  (144).  Hence,  by  applying  the  triangle  inequality  and  the  infor¬ 
mation  from  (142),  (133),  (152),  and  (135)  consecutively,  we  have 


ll-ffM  - -RpMllr  =  IIPh  -  r«j 

<  lir«  -  rifii  +  iiry  -  rRj 

oo 

<  Yi  IM +  «*,(*)  -  fftMllr 

n=Mk  +  l 


< 


1 

2  k~1 


£j  +  S 


(*) 

P+1 


< 


1 


=  e. 


(153) 


This  establishes  (140). 

Remarks.  Algorithm  III  can  be  easily  adapted  and  modified  to  produce  a  possibly  lower 
degree  rational  model  Rp(z )  that  satisfies  the  design  criterion  (119).  In  fact,  as  a  conse¬ 
quence  of  Theorem  5,  we  can  theoretically  obtain  the  lowest  degree  optimal  Hankel-norm 
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rational  model  Rp(z)  in  the  limiting  case.  One  suggestion  is  to  repeat  this  algorithm  for 
various  values  of  £j,  such  that 


e 

e\  =  — e,  t  =  -  1,  (154) 

L 

for  a  suitably  large  integer  L.  As  to  the  modification  of  the  algorithm  itself,  one  may  wish 
to  try  various  values  of  £,-,  i  =  2, 3, . . . ,  k,  instead  of  £2  =  | £1 ,  £3  =  ^£2  =  jf£i  ,...,£*  = 
^£*_i  =  jrrr£i  as  suggested  in  (124),  (129),  and  (130).  The  only  restriction  is  that 
£1  >  £2  >  £3  >  •  •  •  >  £*•  The  smaller  the  values  of  £/  —  £e-i,  £  =  2 ,...,fc,  chosen, 

the  better  the  chance  is  to  find  the  lowest  degree  rational  model  Rp(z).  Of  course,  more 
computing  time  is  required. 

We  also  remark  that  although  the  Hankel-norm  specification  in  (119)  is  not  as  desirable 
as  the  supremum  (or  uniform)  norm  specification,  it  is  very  close  to  it,  in  view  of  Nehari’s 
theorem  (i.e.  Theorem  2),  since  the  only  thing  that  can  go  wrong  is  an  H°°  (or  analytic) 
additive  factor  which  only  contributes  to  noncausal  information.  In  addition,  the  Hankel- 
norm  specification  is  more  desirable  than  the  L 2  norm  (or  RMS)  specification,  since  it 
follows  from  (31)  that 


llff(z)  -  <  ll#W  -  KpMllr,  (155) 

where  L2  =  L2(\z\  =  1)  and  H(z)  is  assumed  to  be  causal,  in  the  sense  that  H(z)  =  Ha(z) 
or 


tf(z)  =  ]T/inz",  (156) 

n=l 

where  {/*„}  satisfies  (118).  Hence,  as  a  consequence  of  (154),  if  Rp{z )  satisfies  the  design 
criterion  |!f/(z)  —  Rp(z)j|r  <  £,  it  also  satisfies  the  design  criterion  ||H(z)  —  Rp{z)\\l'2  <  £• 
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APPENDIX 


Since  it  is  important  to  be  able  to  relate  a  rational  model 


where  P(z )  and  Q(z)  are  coprime  polynomials  satisfying 

degree(P)  <  degree(Q)  =  M, 
with  the  infinite  Hankel  matrix 


TH  = 


[hi 

h2 

hi 


(157) 


(158) 


(159) 


we  include  in  this  Appendix  a  proof  of  Kronecker’s  theorem  given  as  Theorem  1  on  page 
2.  To  be  consistent  with  the  notation  in  (6),  let 


7«  —  {hi,  hi+i , . . .) 


denote  the  ith 


column  vector  of  Th  •  We  also  need  the  notation 


hi 


hi+M-l 


hi+M-i 

h,+2M~2 


for  the  M-dimensional  cofactors  of  Th  with  leading  entry  h,.  Then 


(160) 


(161) 


H\i  =  Hm  (162) 

is  the  principle  cofactor  of  T h  of  dimension  M  introduced  in  (7)  and  (91).  We  have  the 
following  preliminary  result. 


Lemma  1.  The  infinite  Hankel  matrix  Th  has  finite  rank  =  k  if  and  only  if  the  first 
k  column  vectors  k  of  T h  are  linearly  independent  and  there  exist  k  numbers 

Cj, . . .  ,Ck  such  that 

k 

7 k+e  =  ^  ^ Ci7«'-M— i,  £  =  1,2,...  .  (163) 

1=1 

Furthermore,  if  rank  Th  —  k  <  oo,  then  the  principal  cofactor  H\  =  Hk  is  a  nonsingular 
square  matrix. 

Proof.  By  definition,  rank  Th  <  oo  if  and  only  if  Tp  has  only  a  finite  number  of  linearly 
independent  columns.  Suppose  rank  T h  <  oo.  Let  7i,  - . .  ,7r  be  linearly  dependent;  and 
let  k  be  the  largest  integer,  1  <  k  <  r,  such  that  7i , . . . ,  7*  are  linearly  independent.  Then 
since  71 , . . . ,  7*+i  are  linearly  dependent,  there  exist  constants  cj , . . . ,  c*  such  that 

k 

7Jfc+i=^c,7,.  (164) 

i=  1 
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This  is  (163)  for  £  =  1.  Note  that  both  sides  of  (164)  are  infinite  dimensional  vectors.  By 
deleting  the  first  (£  —  1)  entries  of  these  vectors,  (164)  becomes  (163)  for  £  =  2,3,...  . 
That  is,  each  tth  column  of  T//,  for  i  >  k,  is  a  linear  combination  of  its  previous  k 
columns;  and  hence,  by  using  (163)  repeatedly,  each  7i,  i  >  k,  is  a  linear  combination  of 
7i,  •  •  • ,  7*.  So,  for  finite  rank  Th,  its  rank  is  given  by  the  largest  k  for  which  71 , . . . ,  7* 
axe  linearly  independent.  Of  course,  rank  Th  =  00  if  and  only  if  such  a  k  does  not  exist. 
This  proves  the  first  statement  in  the  lemma.  Now,  suppose  that  rank  Th  —  k  <  00. 
By  repeated  applications  of  (163),  it  is  clear  that  every  minor  det  Hlk,  i  =  2, 3, . . . ,  is  a 
constant  multiple  of  the  principle  minor  det  H\.  Recall  that  if  rank  Th  =  k,  then  F h  has 
some  k  dimensional  cofactor  with  nonzero  determinant,  (i.e.  det  H'k  =£  0  for  some  i),  so 
that  det  Hk  =  det  Hi  ^  0. 

We  are  now  ready  to  prove  Theorem  1.  Suppose  that  H(z)  in  (157)  satisfies  (158). 
Then  by  dividing  both  the  numerator  and  denominator  by  the  leading  coefficient  of  Q{z), 
we  may  write 


{ 


P(z)  =  b1zM~i  +  ---  +  bM 
Q{z)  =  z^1  +  a\zM  1  +  ■  •  •  +  gjw 


61  =  hi 

^2  —  ^2  4*  hjOj 


6jW  —  /lAf  +  hfri -lai  +  •  •  •  +  hidM-i 

and  by  equating  the  coefficients  of  z-1 ,  z~2, . . .,  we  have 


M 

hM+l  +  =  0 

i=l 

M 

h\1+ 2  +  ]C  =  0 

i=l 


Now,  if  we  define 

Ci  —  gm — »+i  1  *  =  1,  •  •  • ,  JW , 

then  we  observe  that  (168)  is  equivalent  to 


M 


7M+1  =  y  -g.7 M-i+i 

i=l 

M 

=  X>7,  , 


(165) 


so  that  (157)  is  equivalent  to 

( h\z  1  +  /122  2  +  ■  •  ')(2^  4"  dyz^1  1  +  •  •  •  4~  ojvf )  =  5] zM  1  +  •  •  •  +  bm •  (166) 

Hence,  by  equating  the  coefficients  of  zM_1, . . . ,  z,  1,  we  have 


(167) 


(168) 


(169) 


(170) 


1=1 
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which,  by  the  same  argument  as  the  proof  of  Lemma  1,  yields 

M 

lM+t  -  £  =  1,2,....  (171) 

i=i 

Since  P(z )  and  Q(z )  are  coprime  and  the  leading  coefficient  of  Q(z )  is  normalized  to  be 
1,  the  set  of  coefficients  {a;},  and  hence  {cj},  in  (170)  is  unique.  That  is,  71, . . .  ,7m  are 
linearly  independent.  Hence,  by  Lemma  1,  rank  Th  =  M. 

Conversely,  suppose  that  rank  Tn  =  M .  Then  by  Lemma  1,  we  can  find  coefficients 
cii  •  ■  ■  ,cm,  such  that  (171)  holds.  Hence,  defining  ai, . . .  ,a\i  and  6j, ...  ,6m  by  (169)  and 
(167),  respectively,  we  have  both  (167)  and  (168),  which  yields  (166);  or  equivalently,  we 
have  H(z )  =  P(z)/Q(z)  by  using  (165)  to  define  P(z)  and  Q(z).  The  linear  independence 
of  7i,...,7 m  is  equivalent  to  the  uniqueness  of  {ai)  and  {6i},  i  —  1  ,...,M,  which,  in 
turn,  implies  that  P(z)  and  Q(z)  are  coprime. 
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